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Summary 


In  the  last  year  the  work  on  this  grant  has  involved  the  efforts  of  Eric  Sheppard  in  pursuit  of  his 
doctoral  degree.  Additionally  Jean-Marc  Chanty  was  supported  for  a  time  on  this  grant.  The  work 
of  Eric  Sheppard  is  summarized  in  the  following  attachment.  The  work  of  Jean-Marc  Chanty  is  a 
continuation  of  the  research  presented  in  the  attached  paper.  In  addition  Eli  Niewood  received  partial 
support  from  this  grant.  His  work  has  resulted  in  a  paper  to  be  submitted  for  publication.  A  summary 
of  his  work  and  the  paper  is  attached. 


2 


Nonequilibrium  and  Radiation  in  MPD  Plasmas 


AFOSR  Grant  86-0019 
Manuel  Martinez-Sanchez 

Dept  of  Aeronautics  and  Astronautics 
MIT, Cambridge,  MA  02139 


This  is  a  summary  of  work  contributed  by  Eric  J.  Sheppard,  graduate  student,  under  the 
AFOSR  grant.  Some  early  results  were  reported  in  “Spectroscopic  Investigation  of  the  Exit  Plane 
of  an  MPD  Thruster”  (Kilfoyle,  Martinez,  Heimerdinger  and  Sheppard;  IEPC  paper  #  88-027). 

The  basic  co-axial  MPD  thruster  is  shown  in  figure  1.  The  plasma  is  accelerated  primarily  by 
the  Lorenz  force,  and  is  heated  by  Ohmic  dissipation.  For  this  analysis,  a  1-D  approach  is  taken, 
represented  by  the  idealized  set-up  shown  in  figure  2.  That  is,  in  this  work,  the  flow  through  a 
long,  thin  channel  is  investigated,  with  particular  interest  in  the  radiative  emission  from  the  plasma. 
Figures  3  and  4  show  some  of  the  detailed  physics  that  have  been  included  in  models  at  MIT;  most 
of  these  will  be  in  the  final  model  developed  by  this  work.  The  areas  of  focus  here  are  radiation, 
inelastic  collisions,  and  ambipolar  diffusion.  Following  is  a  discussion  of  work  done,  and  under  way. 

A  finite-rate  approcoh  to  the  kinetics  con  .'rn’  ig  level  populations  in  atoms  was  motivated  by 
the  earlier  experimental  work  of  Heimerdinger  ^  .  J.  Heimerdinger,  Fluid  Mechanics  in  a  Mag- 
netoplasmadynamic  Thruster,  Doctoral  dissertation,  Massachusetts  Institute  of  Technology,  Dept, 
of  Aeronautics  and  Astronautics,  1988)  and  Kilfoyle  (D.  B.  Kilfoyle,  Spectroscopic  Analysis  of  a 
Magnetoplasmadynamic  Arcjet,  Masters  thesis,  M.I.T.,  Dept,  of  Aero.  &  Astro.,  1988),  who  took 
spectroscopic  measurements  at  the  exit  plane  of  an  MPD  thruster.  There  is  interest  both  in  the 
detailed  population  rates  and  the  spectral  intensities  of  line  radiation  emitted  from  electronically 
excited  atoms. 

First,  the  effect  of  inelastic  collisions  (excitation/de-excitation,  ionization/recombination),  line 
radiation,  and  ambipolar  diffusion  (to  the  walls,  where  recombination  to  the  ground  occurs)  on  the 
populations  of  the  excited  energy  levels  in  both  argon  (atom  and  ion)  and  hydrogen  was  investigated. 
This  is  referred  to  as  the  CRD  model,  and  is  also  used  for  calculating  ionization  fractions.  The 
initial  CRD  model  is  applicable  to  a  steady-state,  static  discharge  in  a  plasma,  and  is  based  on 
the  formulation  presented  in  several  sources,  in  particular  Mitchner  and  Kruger  ( Partially  Ionized 
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Gases,  Chapter  9).  Continuity  equations  for  number  densities  n,  with  e  representing  the  ion  (since 
ne  =  rij),  k  =  1  the  ground  state  atom,  and  k  >  1  the  excited  states  are: 


dne 

~dt 


+  V  ■  (netf) 


dnk= i 
dt 


+  V • (nfc=1u) 


gnt>i 

dt 


+  V ‘ (rti>iu) 


L2 

L2 


nk>i 


For  the  steady-state,  static  case,  the  left  hand  side  of  each  of  these  equations  is  identically  equal 
to  zero.  Da  is  the  ambipolar  diffusion  coefficient,  L  is  an  equivalent  diffusion  radius,  and  the  rate 
equations  are: 


—  Me  )  ’  nk$kc  ^  d"  Ackfick) 

k  k 

nk  =  ne  ^2njSjk  ~  nfc(5Z(««5fcy  +  AkjPkj)  +  +  Skc)} 

]<k  j<k  l>k 

+  ^i{neSik  +  Atkfiik)  +  n] ( neSck  +  Ack0ck) 

l>k 

Here,  the  5  terms  are  calculated  collisions!  rate  terms  and  the  A  terms  are  transition  probabilities 
for  the  radiation  events.  The  ratio  of  radiation  lost  from  the  volume  to  the  total  radiation  emitted  is 
denoted  by  the  radiative  escape  factor,  fi,  which  ranges  from  0  (optically  thick)  to  1  (optically  thin). 
In  all  of  the  terms,  the  subscripts  refer  to  processes  which  involve  a  change  from  the  level  denoted 
by  the  first  index  to  the  second  (c  here  refers  to  the  continuum,  or  ion  state).  It  is  important  to 
note  that  this  is  not  a  self-consistent  problem  as  stated  -  the  temperatures  are  reasonably  chosen, 
but  still  arbitrary,  as  is  the  total  number  density.  Also,  only  k  equations  are  solved  for  the  k  +  1 
populations  (the  k  levels  of  the  atom  and  the  ion/electron;  the  ground  state  ( k  =  l)  equation 
is  dropped)  -  the  electron  density  is  arrived  at  iteratively,  as  it  appears  as  a  cubic  term  in  the 
equations. 

Results  show  that  radiation  and  diffusion  effects  on  the  populations  of  excited  levels  are  largest 
(relative  to  a  purely  collisional  model)  at  low  temperatures  and  number  densities,  with  diffusion 
(across  a  typical  2  cm  MPD  channel)  having  the  larger  influence.  Ionization  fractions  may  be  lower 
than  thermal  equilibrium  results  by  several  orders  of  magnitude  under  these  conditions.  Figure  5 
shows  typical  population  results  on  a  Boltzmann  plot  for  argon.  The  solid  line  represents  what 
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the  populations  curve  would  look  like  if  all  the  levels  were  in  Saha  equilibrium.  The  slope  here  is 
an  indicator  of  the  electron  temperature,  and  is  in  general  agreement  with  the  equilibrium  curve. 
However,  for  smaller  diffusion  radii,  the  mid-range  of  energy  levels  may  not  reflect  the  proper 
temerature.  This  difference  is  even  more  dramatic  in  the  argon  atom.  In  both,  only  the  uppermost 
levels  are  always  in  equilibrium  with  the  electrons. 

Next,  a  2-gas  static  model  was  developed.  This  work  is  designed  to  follow  up  on  a  procedure 
used  by  Kilfoyle  and  Heimerdinger.  The  argon  plasma  was  seeded  with  hydrogen  in  an  attempt  to 
use  hydrogen  lines  for  the  spectrographic  diagnostics.  A  2-gas,  2-level  each  model  is  in  use.  The 
formulation  of  the  problem  is  similar  to  the  1-gas  static  model,  with  the  only  interaction  between 
the  two  gases  coming  from  electronic  collisions.  The  individual  n*  formulations  are  the  same  as 
used  earlier,  with  the  complete  set  of  equations: 


n;jl  +  n*‘  +  n 
dnl 


i  k 

,  al 


ntot 


+  V  •  (nf  C) 

r  +  v '(»,*’«) 


dt 

dn’2 

at 

—gf1  +  V  •  (.Si,") 


+  V-«1u)  = 


•at  VW 

n<  -~u~ 

■s2  D?n« 
ne  JT~ 


■si 

nJ>l 


•  s2 
nk>  1 


i1  -  A)(X]R>1  +  fn‘)  =  x(Y,nk  +  i1  -  0ne) 


Here,  si  and  s2  denote  the  two  gases  used.  The  factor  \  is  the  fraction:  £j~,  which  is  a  user 
input,  and  £  is  which  is  solved  for  along  with  the  populations. 

At  some  temperatures  the  ionization  fraction  calculated  for  hydrogen  was  smaller  than  that  for 
argon,  and  at  the  number  densities  looked  at  (around  1015cm~3)  the  difference  between  the  two 
was  not  very  large.  This  seemed  counter-intuitive  at  first,  but  it  turned  out  to  be  due  to  what 
is  referred  to  here  as  the  crossover  temperature.  That  is,  the  difference  in  the  degeneracies  of  H 
and  A  mean  that  there  is  a  temperature  (around  10100  K)  where,  for  any  njo(,  the  (Saha  thermal) 
equilibrium  ionization  fractions  are  equal:  =  a//.  Below  this  crossover  temperature,  a a  < 

and  above  it,  >  a//.  The  crossover  temperature,  Tx,  can  be  found  from  the  Saha  equation  by 
setting  aSl  =  a,a  for  any  two  gases  si  and  S2  (where  9  =  9 's  a  degeneracy,  »  represents  the 

ion,  and  1  the  ground  state  atom): 
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%F^-  =  lr>([-n-n 

*  x  9 1  9i 

Figure  6  shows  the  ionization  fractions  of  A  and  H  as  functions  of  the  electron  temperature,  calcu¬ 
lated  by  the  Saha  equilibrium  equation.  The  vertical  line  is  drawn  at  about  the  theoretical  crossover 
temperature  for  A  and  H,  and  the  total  number  density  increases  by  one  order  of  magnitude  with 
every  pair  of  lines  moving  to  the  right. 

Results  from  the  2-gas  model  for  low  mass  fractions  (hydrogen:total)  show  that  the  hydrogen 
acts  as  a  fairly  poor  seed  (it  doesn’t  increase  the  ionization  fraction  much)  at  the  temperatures 
below  Tx,  and  in  fact  it  lowers  a  a  little  at  above  Tx.  This  may  be  good,  assuming  that  the  only 
interaction  between  the  two  gases  is  between  their  electrons,  because  the  addition  of  the  hydrogen 
has  a  small  effect  on  the  bulk  plasma  properties  (i.e.,  conductivity).  However,  we  have  the  possible 
advantage  of  gaining  access  to  hydrogen  lines  for  the  spectrometer. 

Channel-flow  physics  need  to  be  included  to  study  the  effects  of  a  plasma  out  of  thermal  equi¬ 
librium  moving  through  a  discharge  (i.e.,  what  we  see  in  an  MPD  thruster).  This  makes  the  work 
self-consistent;  that  is,  the  temperatures  and  total  density  are  no  longer  arbitrarily  input,  but 
are  calculated  in  the  channel.  A  1-D  flow  model  (a  space-marching  scheme  following  Niewood’s 
work  (E.  Niewood,  Transient  One  Dimensional  Numerical  Simulation  of  Magnetoplasmadynamic 
Thrusters,  Masters  thesis,  MIT,  Dept,  of  Aero.  &  Astro.,  1989))  is  being  developed.  This  is  a 
two-temperature  (one  temperature  for  the  light  electrons  (Te),  one  temperature  for  the  heavies: 
the  ions  and  neutrals  (Tg)),  one  speed  model,  which  includes  ionization  and  recombination  and 
area  variation.  Radiative  decay,  and  detailed  continuity  equations  for  the  excited  levels  have  been 
added.  The  steady-state  1-D  equations  in  non-dimensional  form  are  (£  is  the  nondimensional  axial 
coordinate)  : 


^ua  =  1 


dot  , 

—  =  a.(a  -  D) 


dpi 

d£ 


=  a(P  +  D) 


d0k>  l 

d£ 


a0k>\ 
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du  d  r  .2i 
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daTe  dp  2 J2  . 

P—77-  -  (if  -  l)aT'17  =  (l~  +  — -j 


d£ 


dZ 


RmU1 


dT0  .  dp  ,  R 

~  [1  ~  = (1  “  '>— 

Where  r  is  the  viscous  force;  $  is  the  viscous  dissipation;  a  is  the  ionization  fraction  (^p *-) ; 
D  is  the  ambipolar  diffusion;  /?*  =  ^7,  where  n*  is  the  number  density  of  the  level  (level  1  is 
the  ground  state);  R  is  the  radiative  energy  loss  from  excited  species;  and  Rm  =  pt0auL  is  the 
magnetic  Reynolds  number.  The  a-like  terms  represent  the  net  kinetic  rate  of  population  of  the 
level  involved,  and  it  is  assumed  that  diffusion  to  the  walls  results  in  a  recombination  event  resulting 
in  a  ground  state  atom.  Using  these  equations,  along  with  the  definition  of  the  speed  of  sound  (c„) 
and  Mach  number  (Af),  an  equation  of  state,  and  Ohm’s  law,  the  following  expression  for  the 
velocity  gradient  is  found: 


i«  +  (T  -  1)(*^  -  +  £("*- 

d£  A/2  -  1  1  j 

This  form  poses  a  numerical  difficulty  when  integrating  from  the  inlet,  in  that  both  the  numera¬ 
tor  and  denominator  must  go  to  zero  simultaneously.  This  problem  is  avoided  by  using  L’Hopital’s 
rule  to  evaluate  the  derivative  at  the  sonic  point  (Af  =  1),  and  integrating  forward  to  the  exit 
plane.  If  the  proper  condition  is  met  there,  namely,  6  =  0,  then  a  back-integration  is  initiated  from 
the  sonic  point  to  obtain  the  inlet  conditions.  At  present,  this  is  being  carried  out  for  a  2-level 
(one  excited  state  and  the  ground)  argon  atom;  some  results  are  shown  in  Figure  7.  This  same 
approach  was  used  earlier  in  a  one-fluid  (MHD)  model,  and  figure  8  shows  the  variation  of  the 
sonic  passage  point  with  the  magnetic  Reynolds  number.  This  behavior  validates  the  inner-outer 
expansion  used  by  Martinez  in  an  earlier  publication.  (The  Structure  of  Self-Field  Accelerated 
Flows”,  AIAA-87-1065) 

The  expediency  of  this  1-D  approach  has  a  shortcoming;  all  the  transverse  gradients,  which  are 
critical  for  radiative  transfer,  must  be  approximated  algebraically.  To  this  end,  a  brief  examination 
of  boundary  layers  in  the  MPD  will  be  made.  Analytical  work  has  started,  and  a  simple  numerical 
model  is  being  tested. 
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The  final  step  is  to  analyze  the  spectral  intensities  of  hydrogen  and  argon  lines;  this  to  be  com¬ 
pared  to  previous  experimental  work,  and  used  in  guiding  future  experimental  diagnostic  methods. 
In  particular,  the  choices  of  seeding  options  and  spectral  lines  to  observe,  and  the  relationship  be¬ 
tween  the  intensities  measured  and  the  temperatures  and  (derived)  number  densities  of  the  plasma 
will  be  studied.  Thus,  the  result  of  this  study  will  be  a  1-D  numerical  model  (axial  flow  with 
transverse  effects  either  assumed  or  approximated  in  algebraic  form)  for  an  argon  plasma  seeded 
with  hydrogen  in  an  MPD  thruster.  Multi-level  models  will  be  used  in  order  to  capture  the  details 
of  the  radiating  excited  levels  of  interest. 

Preliminary  work  agrees  with  the  results  of  Niewood,  that  there  are  two  distinct  temperatures, 
and  that  the  ionization  rate  terms  do  have  a  significant  effect  on  the  flow  variables.  Future  work 
will  extend  the  existing  model,  as  well  as  add  in  more  of  the  processes  shown  in  figures  3  and  4. 
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Two  Fluid  Numerical  Simulations  of  MPD  Plasmas 
AFOSR  GRANT-86-0019 
Manuel  Martinez-Sanchez 

Dept  of  Aeronautics  and  Astronautics 
MIT,Cambridge,  MA  02139 


This  is  a  summary  of  work  done  by  Eliahu  Niewood,  a  graduate  student,  supported  in  part  by  the  AFOSR 
grant. 


Technical  Discussion 

A  variety  of  different  types  of  electric  propulsion  have  been  put  forward  as  appropriate  for  space 
applications.  One  class  of  electric  propulsion  device  is  the  magnetoplasmadynamic,  or  MPD,  thruster. 
MPD  thrusters  use  the  Lorentz  force  produced  by  charged  particles  moving  in  an  electric  field  to  accelerate 
the  propulsive  fluid  and  produce  thrust. 

MPD  channels,  although  simple  to  construct,  are  hard  to  analyze  because  of  the  complex  interaction 
between  the  magnetic  field  and  the  fluid.  One  way  to  model  these  devices  and  to  predict  their  performance  is 
to  simulate  them  numerically.  Some  numerical  work  has  already  been  done  with  MPD  thrusters.  Martinez 

[5]  and  Kuriki  [3]  study  a  one  fluid,  fully  ionized  model  of  a  steady  quasi  one  dimensional  flow.  Minakuchi 

[6]  includes  the  effect  of  axial  heat  conduction  in  his  steady,  one  fluid,  quasi  one  dimensional  model. 
Subramaniam  [7]  and  Lawless  [4]  worked  with  both  the  one  fluid  model  and  a  partly  ionized  two  fluid 
model.  Heimerdinger  [lj  also  works  with  a  partly  ionized,  thermal  equilibrium  model,  which  includes 
viscosity.  Other  work  has  been  done  with  two  dimensional  models. 

This  research  uses  a  quasi  one  dimensional  two  fluid  model.  The  model  consists  of  fluid  equations 
for  the  total  density,  the  ionization  fraction,  the  axial  velocity,  the  electron  temperature,  and  the  heavy 
soecies  temperature,  as  well  as  a  magnetic  field  equation  derived  by  combining  Maxwell’s  equations  with 
ohm's  law.  It  is  assumed  that  all  the  particles  have  the  same  axial  velocity,  and  that  the  ion  and  neutral 
temperatures  are  the  same. 

A  number  of  source  terms  are  used  to  model  various  phenomena  in  the  flow.  The  electron  temperature 
equation  includes  axial  heat  conduction.  A  number  of  source  terms  are  used  to  model  various  phenomena 
in  the  flow.  The  Hinnov-Hershberg  model  is  used  to  determine  the  ionization  and  recombination  at  each 
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axial  location.  Ambipolar  diffusion  is  included  by  assuming  a  parabolic  number  density  distribution  across 
the  channel.  The  contribution  of  viscosity  to  the  momentum  and  the  heavy  species  temperature  equations 
is  also  included  by  assuming  a  parabolic  velocity  distribution.  Collisional  energy  transfer  between  electrons 
and  heavy  species  is  also  modeled.  The  electrical  conductivity  at  each  axial  location  is  computed  using 
the  Spitzer-Harm  formula.  The  viscosity  and  heat  conduction  coefficients  are  also  computed  at  each  axial 
location.  All  of  the  equations  are  written  in  their  unstead}  forms,  and  the  numerical  method  is  an  unsteady 
method,  although  only  steady  state  results  have  been  examined  so  far. 

The  numerical  method  used  is  a  combination  of  various  schemes.  This  is  due  to  the  different  types  of 
equations  and  their  different  time  scales.  The  magnetic  field  equation  and  the  heat  conduction  term  in 
the  electron  temperature  equation  are  integrated  using  McCormack’s  method  a  number  of  times  for  each 
step  of  the  other  equations.  The  remaining  equations  are  stepped  using  either  the  Donor  Cell  method  or 
a  modified  Rusanov  method. 

A  number  of  the  inlet  boundary  conditions  of  the  flow  are  determined  by  the  physical  characteristics  of 
the  flow,  the  mass  flow  per  unit  area,  the  applied  current,  and  the  total  inlet  enthalpy.  The  inlet  ionization 
fraction  is  physically  zero,  but,  because  of  the  limitations  of  the  model,  the  inlet  ionization  fraction  is 
set  to  some  small  number.  The  inlet  heavy  species  temperature  is  set  to  room  temperature.  For  cases 
with  heat  conduction,  the  electron  temperature  gradient  at  the  inlet  is  set  to  zero.  For  cases  without  heat 
conduction,  the  inlet  density  is  computed  using  a  downwind  difference. 

The  exit  boundary  conditions  for  supersonic  flow  are  zero  gradient  in  all  variables  except  the  magnetic 
field,  which  is  set  to  zero  at  the  exit,  and  the  heavy  species  temperature,  which  is  chosen  so  that  there  is 
zero  gradient  in  the  sum  of  the  fluid  and  magnetic  pressures.  For  subsonic  flow  the  pressure  is  set  to  some 
small  external  pressure,  and  the  other  variable  are  chosen  so  that  there  is  zero  gradient  in  the  Riemann 
invariants. 

The  method  was  tested  by  simplifying  the  equations  to  be  the  unsteady  analogue  of  the  equations 
used  by  Martinez  [5],  Running  the  simulation  to  steady  state  provided  results  similar  to  those  found  by 
Martinez  for  a  number  of  different  magnetic  Reynolds  numbers  and  for  both  constant  and  variable  area 
channels. 

The  full  set  of  equations  were  then  used  to  examine  the  effect  of  various  phenomena  on  the  flow  and 
performance  characteristics  of  the  thruster.  The  thruster  was  taken  to  be  20  cm  long  with  an  electrode 
seperation  of  2  cm.  The  inlet  ionization  and  total  enthalpy  were  chosen  to  be  small  enough  so  that  their 
effect  on  the  flow  was  negligible. 

First  ambipolar  diffusion  was  added  to  the  baseline  case,  then  viscosity,  and  finally,  heat  conduction. 
As  shown  in  Table  1,  ambipolar  diffusion  and  heat  conduction  had  little  effect  on  performance.  Viscosity 
however  played  a  major  role  in  determining  thruster  performance  and  flow  characteristics.  Viscosity  causes 
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Case 

Bmag 
or  Bo 

ambipolar 

diffusion 

viscosity 

heat 

conduction 

channel 

E, 

V  Jm 

T 

N/m 2 

■ 

1-fluid 

20 

no 

no 

no 

CAC 

345 

4370 

0.69 

1-fluid 

4.9258 

no 

no 

no 

CAC 

433 

4630 

0.62 

2-fluid 

0.1 

no 

no 

no 

CAC 

399 

4256 

0.57 

2-fluid 

0.1 

yes 

no 

no 

CAC 

401 

4180 

0.55 

2-fluid 

0.1 

yes 

yes 

no 

CAC 

258 

2100 

0.21 

2-fluid 

0.1 

yes 

yes 

yes 

CAC 

287 

2300 

0.23 

2-fluid 

0.15 

yes 

yes 

yes 

CAC 

519 

4550 

0.33 

2-fluid 

0.2 

yes 

yes 

yes 

CAC 

860 

6980 

0.33 

2-fluid 

0.1 

yes 

yes 

FFC 

340 

3053 

0.34 

Table  1:  Thrust  and  Efficiency  for  All  Cases 


the  heavy  species  temperature  to  increase  rapidly  to  levels  of  1  to  10  eV  and  causes  the  flow  to  be  frictionally 
choked  with  exit  velocities  on  the  order  of  5000  meters/second  . 

It  was  also  desired  to  determine  the  effect  of  increasing  the  total  current  supplied  to  the  thruster, 
and  hence  the  inlet  magnetic  field.  Increasing  the  total  current  from  80  kAmp/meter-depth  to  120  kA/m 
caused  the  flow  to  develop  a  large  current  concentration  at  the  exit,  with  a  corresponding  increase  in  the  exit 
electron  temperature  and  ionization.  Increasing  the  total  current  to  160  kA/m  caused  this  concentration 
to  become  so  large  that  the  simulation  did  not  converge  to  a  steady  state.  As  shown  in  Table  1,  increasing 
the  total  current  to  120  kA/m  resulted  in  an  increase  in  the  thruster  efficiency. 

It  is  possible  that  the  large  exit  current  concentrations  in  the  higher  total  current  cases  is  due  in  part 
to  the  static  instability  described  by  Heimerdinger  [2].  Physically  this  instability  would  first  appear  near 
the  inlet  and  then  travel  down  the  channel.  When  the  current  buildup  reaches  the  exit,  the  current  would 
bend  out  of  the  channel  and  the  concentration  would  dissipate  .  The  current  concetration  would  then 
reappear  at  the  beginning  of  the  channel.  However,  in  a  one  dimensional  model  with  zero  exit  current,  the 
instability  is  forced  to  remain  at  the  exit. 

The  last  effect  examined  was  the  effect  of  area  variation.  The  results  for  the  full  model  at  the  80  kA/m 
current  for  a  constant  area  channel  were  compared  to  the  results  for  a  parabolic  channel  with  an  throat 
area  of  2  cm  ,  an  inlet  to  throat  area  ratio  of  1.5,  and  an  exit  to  throat  area  ratio  of  2.  Area  variation 
alleviated  the  exit  current  buildup  somewhat  and  resulted  in  an  increase  in  performance. 

As  detailed  earlier,  a  model  with  a  one  dimensional  magnetic  field  which  goes  to  zero  at  the  channel 
exit  is  problematic.  To  alleviate  these  problems,  a  two  dimensional  magnetic  field  simulation,  coupled  with 
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a  quasi  one  dimensional  fluid  calculation,  was  developed.  This  simulation  allows  the  current  to  bend  at 
the  exit.  However,  the  quenching  of  the  exit  current  has  not  yet  been  simulated  numerically. 

Some  work  has  also  been  done  in  an  attempt  to  combine  the  one  dimensional  MPD  simulation  with  a 
model  of  anomalous  dissipation  in  MPD  thrusters.  The  dissipation  arises  due  to  the  presence  of  a  modified 
two  stream  instability  in  the  plasma.  The  one  dimensional  code  produces  the  parameters  at  each  axial 
location  needed  by  the  dissipation  model  to  determine  if  the  instability  exists  and  what  its  effect  will  be,  in 
terms  of  a  modified  conductivity  and  electron  and  ion  heating  rates.  The  one  dimensional  code  will  then 
take  into  account  the  effects  of  the  instability  .  By  iterating  this  process  a  number  of  times,  a  steady  state 
should  be  reached. 
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Abstract 

Magnetoplasmadynamic,  or  MPD,  thrusters  are  a  promising  method  of  propulsion  for  a  variety 
of  different  space  missions.  This  research  develops  and  analyzes  a  numerical  simulation  of  a  quasi 
one  dimensional  model  for  an  MPD  thruster.  A  finite  difference  scheme  is  used  to  integrate  the  fluid 
equations  for  each  species  and  a  magnetic  field  equation  derived  from  Maxwell’s  laws.  The  model 
includes  separate  electron  and  heavy  species  temperatures,  varying  conductivity,  varying  ionization 
fraction,  collisional  energy  transfer  between  heavy  particles  and  electrons,  averaged  viscosity  and 
ambipolar  diffusion,  and  electron  heat  conduction.  Both  constant  area  and  variable  area  channels 
are  examined.  The  applied  current  in  the  cases  studied  ranges  from  79.6  —  to  159  mV^r7fp,T 

for  an  inlet  mass  flow  of  0.5  ■  •  It  is  shown  that  thermal  equilibrium  is  not  a  valid  assumption 

in  a  typical  MPD  thruster.  It  is  also  found  that  viscosity  plays  a  significant  role  in  determining 
thruster  performance.  Area  variation  is  also  found  to  have  a  significant  effect  on  performance. 


Nomenclature 

Symbols 

A  Channel  Area 
a  Ionisation  fraction 
B  Magnetic  field  strength 
Da  Ambipolar  diffusion  coefficient 

Ei  Elastic  momentum  transfer  between  electrons  and  heavy  species 
E,  Ionization  energy 
<  Electric  charge  of  a  proton 
c0  Permittivity  of  vacuum 
T  Spitzer  logarithm 
H  Channel  height 
/  Total  electric  current 
J  Electric  current  density 
K  Heat  conduction  coefficient 
H  Viscosity  coefficient 

Permeability  of  vacuum 
n  Number  density 
n«  Ionization  rate 

v„  Collision  frequency  of  species  s  with  species  r 
P  Scalar  fluid  pressure 

Q,r  Collision  cross  section  of  species  s  with  species  r 
p  Mass  density 
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51  Source  term  for  continuity  equation 

5 2  Source  term  for  momentum  equation 

53  Source  term  for  energy  equation 
a  Scalar  electrical  conductivity 

T  Temperature 
U  Average  velocity  of  fluid 

Subscripts 

a9=i,e,n  Value  of  a  for  ion,  electron,  or  neutral  species 


1  Introduction 

A  variety  of  different  types  of  electric  propulsion  have  been  put  forward  as  appropriate  for  space 
applications.  One  class  of  electric  propulsion  device  is  the  magnetoplasmadynamic,  or  MPD  thruster, 
which  utilizes  the  Lorentz  force  produced  by  charged  particles  moving  in  a  magnetic  field  to  accelerate 
the  propulsive  fluid  and  produce  thrust. 

MPD  channels,  because  of  the  interaction  between  the  magnetic  field  and  the  fluid,  are  extremely 
complicated  and  hard  to  analyze.  One  possible  way  to  understand  these  devices,  and  to  predict  their 
performance,  is  to  simulate  them  numerically.  Some  numerical  work  has  already  been  done  with  MPD 
thrusters.  A  number  of  works  deal  with  steady  quasi  one  dimensional  flow.  Martinez  [14]  and  Kuriki 
[12]  study  a  one  fluid,  fully  ionized  model.  Minakuchi  [17]  includes  the  effect  of  heat  conduction. 
Subramaniam  [20]  and  Lawless  [13]  worked  with  both  the  one  fluid  model  and  a  partly  ionized  two  fluid 
model  with  thermal  equilibrium.  They  include  both  viscosity  and  heat  conduction  to  the  wall  in  their 
model.  Heimerdinger  [7]  also  works  with  a  partly  ionized,  thermal  equilibrium  model,  which  includes 
viscosity.  Auweter-Kurtz  [l]  works  with  a  quasi  one  dimensional  model  which  allows  for  seperate  heavy 
and  electron  temperatures,  and  varying  ionization  fraction,  in  parts  of  the  thruster.  Other  existing 
works,  by  Chanty  [2],  Sleziona  [  19] ,  and  Park  [18]  examine  fully  two  dimensional  flow  in  both  steady 
and  unsteady  cases,  again  with  a  one  fluid,  fully  ionized  model. 

The  model  used  in  this  research  is  among  the  first  to  study  two  fluid  flow  throughout  the  thruster 
with  separate  heavy  and  electron  temperatures.  It  also  includes  a  number  of  effects,  such  as  ambipolar 
diffusion  and  collisional  energy  transfer,  not  studied  previously. 

2  Model 

2.1  Assumptions 

Due  to  the  complexity  of  the  problem,  the  model  includes  a  number  of  simplifying  assumptions. 
These  are 

1.  The  flow  is  quasi  one  dimensional. 

2.  Ions  and  neutrals  are  tightly  coupled,  so  that  they  have  equal  temperatures. 

3.  All  particles  have  the  same  average  velocity  in  the  axial  direction. 

4.  The  plasma  is  quasi  neutral  with  singly  ionized  ions  only. 

5.  No  axial  current  (the  Hall  effect  is  ignored). 

The  working  fluid  in  this  model  is  Argon,  with  a  molecular  mass  of  6.525  x  10-26  and  an 

ionization  energy  of  2.53  x  10-18  ■ 
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2.2  Governing  Equations 

For  a  channel  wirh  varying  area,  where  A  =  HD  =  H-  (unit  depth),  the  governing  equations  are 


Global  Continuity: 

dp  A  dpUA 
dt  dz 

(i) 

Global  Momentum: 

dpUA 

dt 

+  a*?A+A>{p+&).M 
dz  dz 

(2) 

Electron  Density: 

dpocA  dpaUA 
a.  +  a  =Am,5l 

dt  dz 

(3) 

Electron  Energy: 

dpaTeA 

dt 

dpotUTeA  2  dU A  ,2  rru 

3,  +  yaT'-  a.  =  A,  k  S3‘ 

(4) 

Heavy  Species  Energy: 

dpT,A 

dt 

dpUTgA  2  ^dUA  ,2m,-,, 

«  4"  pig  «  —  A 

dz  3  dz  3  1c 

(5) 

Magnetic  Field: 

dBA  dUBA  1  A  dcr  dA  dB  A  d2B 

dt  1  dz  1  <7 [Aq  <7  dz  dz  dz  crpto  dz 2 

(6) 

where,  assuming  Coulomb  collisions  are  dominant,  a  is  given  by  the  Spitzer-Harm  formula,  a  — 
in  Si,  and  Te  =  1.24  x  lO7^ .  Also,  a  = 


2.3  Source  Terms 


The  source  terms  in  the  electron  density  equation  represent  loss  and  creation  of  electrons.  One 
process  which  contributes  to  this  source  term,  is  the  rate  of  ionization  due  to  electron  impact,  denoted 
by  ne.  This  process  is  evaluated  using  the  Hinnov-Hirshberg  model  of  ionization  [15],  which  gives 


he  =  i?ne(5n„  -  n2] 


where, 


and 


„  1.09  x  10-20  , .  m6 

R  =  -  -  (tn - ) 

ya  Sec 


S  =  2.9  x  1022T^e  (inm  3) 


(7) 

(8) 

(9) 


Another  process  which  contributes  to  the  loss  of  electrons  is  ambipolar  diffusion  to  the  walls  '.’r'-ipolar 
diffusion  is  given  by 

dne  d2ne  12  Dane 
dt  ~  a  dx2  ~  H 2 

where  Da  is  the  ambipolar  diffusivity, 


1 

Qin  (r*n  "h  rie) 


(11) 
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and  where  a  parabolic  distribution  has  been  assumed  for  the  electrons.  Combining  both  terms  yields, 


51  =  ne  — 


12  Dane 

H 2 


(12) 


The  source  term  in  the  momentum  equation  comes  from  the  viscous  forces  exerted  on  the  moving 
fluid  by  the  walls,  so  that 

'atn 

dx  ) 

■  /  W 


H 

where  n,  the  coefficient  of  viscosity  is  given  by  jl6j, 


(13) 


u  = 

1  - 

a 

a 

kg 

.  —  a  + 

'  Ml  .  .  r) 

(1  -  a)“  +  a 

m  •  sec 

) 

mtC7» 

*<j  _  \ 

0. 40G|  4  rra,  )7  ,/^7|  kT,  )  i 

(jn 

kg  \ 

2Qnn 

( Ifl 

m  sec  > 

e*  In  rg 

iin 

m  sec  > 

Qnn 

=  1.7  x  io-18rtf_  ♦ 

(in 

m2) 

Qxn 

=  1.4  x  10~18 

(in 

m2) 

Qa 

c,= 

U  kTq 

V  1 rm. 

(in 

t) 

r„  = 

1.24  x  10 

e*  111  r, 

32n/ikiT} 


(14) 


(in  m2) 


A  parabolic  distribution  is  assumed  for  the  velocity  across  the  channel  and  ^  is  evaluated  at  the  wall 
to  give 

\2Un 


52  =  -- 


H2 


(15) 


One  source  of  enery  in  the  electron  temperature  equation  is  the  Joule  heating,  given  by  Another 
source  of  energy  is  collisional  energy  transfer  between  electrons  and  the  heavy  particles,  denoted  by  Ei, 
where 


Ei  =  5.67  x  10~28(Te  -  Tg)neu, 


watts  , 


(16) 


where, 

(*n  7)  Q'i  =  (mm2)  Ce=  6214.0^  (in  f£) 

Since  the  internal  energy  as  defined  does  not  include  the  ionization  energy  of  the  electrons,  energy  is 
also  lost  when  an  electron  ionizes,  so  there  is  a  loss  equal  to  Eihc.  Also  included  in  this  source  term  is 
the  axial  heat  conduction,  given  by  ^(Jf,.^^),  where  [15], 


K'  = 


1.7142fc27> 


, .  watt  . 

['n7TK] 


(17) 


Radial  heat  conduction  to  the  electrodes  is  assumed  to  be  small  due  to  the  confining  effect  of  the  sheath, 
although  this  asssumption  needs  to  be  examined  in  more  depth.  Including  all  of  these  processes, 


53e  =  —  -  Ei  -  E,he  +  -f  (A.^p) 
cr  dz  oz 


(18) 


The  source  term  for  the  heavy  particle  temperature  equation  includes  the  energy  lost  to  the  electrons, 
described  earlier,  and  the  heat  produced  by  the  viscous  force  applied  by  the  wall  to  the  fluid.  This  heat 
is  given  by  which  is  averaged  over  the  channel  again  assuming  a  parabolic  velocity  distribution. 

Therefore, 

S3g  =  E,  +  (19) 
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3  Numerical  Method 

3.1  Overall  Scheme 

The  overall  numerical  scheme  is  as  follows: 

1.  Evaluate  all  of  the  source  terms. 

2.  Integrate  the  magnetic  field  equation  the  neccessary  number  of  times,  holding  all  other  variables 
constant.  The  magnetic  field  integration  is  done  using  McCormack’s  method. 

3.  Integrate  the  convective  terms  of  the  fluid  equations  using  a  modification  of  Rusanov’s  method,  for 
the  overall  density  and  momentum  equations  and  the  Donor  Cell  method  for  the  electron  density 
and  species  temperature  equations. 

4.  Add  in  the  contribution  from  the  source  terms. 

5.  Integrate  the  diffusive  part  of  the  electron  temperature  equation,  again  using  McCormack’s  method. 

3.2  Inlet  Boundary  Conditions 

The  mass  flow  per  unit  area  was  chosen  to  be  0.5  in  all  cases.  The  magnetic  field  at  the  inlet 
is  determined  by  the  total  applied  current  /,  assumed  to  be  constant,  Bo  =  The  inlet  ionization 

fraction  is  assumed  to  be  constant  and  small,  equal  to  0.01.  The  inlet  density  is  found  by  a  downwind 
difference,  a  one  sided  finite  difference  formulation  of  the  overall  continuity  equation  at  the  zeroth 
point.  The  total  inlet  enthalpy  is  assumed  to  be  2.1  x  105ijj-  in  the  one  fluid  test  cases,  equivalent  to 

the  enthalpy  of  a  room  temperature  gas,  and  5.5  x  10s ^  in  the  two  fluid  cases,  to  allow  for  the  iniet 
ionization. 

For  cases  without  electron  heat  conduction  the  heavy  species  temperature  at  the  inlet  is  assumed  to 
be  constant  and  equal  to  room  temperature,  300  K.  The  electron  temperature  is  then  found  as  a  function 
of  the  known  variables  and  the  total  inlet  enthalpy.  For  cases  with  electron  heat  conduction,  the  electron 
temperature  is  assumed  to  be  the  same  as  that  at  the  next  inside  point,  and  the  gas  temperature  is  then 
found  as  a  function  of  the  other  known  quantities. 

3.3  Outlet  Boundary  Conditions 

At  the  outlet  it  is  assumed  that  if  the  flow  is  supersonic,  the  variables  at  the  exit  are  set  equal 
to  the  variables  at  the  closest  inside  point  with  the  exception  of  the  magnetic  field,  which  is  set  to 
zero.  If  the  flow  is  subsonic,  the  fluid  variables  p  and  U  are  given  by  the  Riemann  invariants,  while  P  is 
determined  by  the  external  pressure,  assumed  to  be  some  very  small  pressure.  The  electron  temperature 
and  ionization  fraction  are  again  taken  from  the  inside  point  and  the  heavy  species  temperature  is  then 
computed  as  in  the  supersonic  case. 

4  Results 

4.1  One  Fluid  Results 

It  is  possible  to  solve  the  one  fluid  equations  with  the  method  developed  for  the  two  fluid  equations. 
The  one  fluid  case  therefore  served  as  a  valuable  test  for  the  method,  because  it  has  already  been 
analyzed  by  others.  Martinez  [14]  uses  a  space  marching  method  with  an  inner-outer  expansion  to 
find  the  steady  state  solution  to  the  one  fluid  equationsand  examines  the  nature  of  the  — *  oo  limit 
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Interelectrode  Separation:  2  cm 

fV  =  0.5  — P— 

Channel  Length  =  20  cm 

Bo  =  0.1  Tesla 

T0  =  300  K 

I  =  79.6  kA.mr 

m  depth 

Table  1:  Thruster  Characteristics 


by  means  of  inner-outer  expansions  at  the  inlet  and  the  exit.  A  test  case  was  run  using  the  thruster 
described  in  Table  1.  The  electric  field  calculated  for  the  test  case  was  within  2%  of  that  calculated  by 
Martinez. 

4.2  Comparison  of  Two  Fluid  and  One  Fluid  Results 

After  the  one  fluid  results  had  been  obtained  and  verified,  the  equivalent  two  fluid  case  was  run  to 
allow  comparisons  of  the  one  fluid  and  two  fluid  models.  This  case  did  not  include  ambipolar  diffusion, 
viscosity  or  heat  conduction,  as  these  effects  are  not  included  in  the  one  fluid  model.  One  important 
feature  of  these  results,  labeled  Case  1  in  Figures  1  -  4,  is  the  large  discrepancy  between  the  electron 
and  heavy  species  temperatures.  The  heavy  species  temperature  is  an  order  of  magnitude  smaller  than 
the  electron  temperature.  This  would  imply  that  the  thermal  equilibrium  assumed  by  the  one  fluid 
model  is  not  a  good  assumption.  This  large  difference  in  temperature  arises  because  there  is  no  effective 
mechanism  in  the  model  of  Case  1  for  heating  the  heavy  species.  The  only  source  term  in  the  heavy 
species  temperature  equation,  collisional  energy  transfer  with  the  electrons,  is  a  relatively  weak  effect. 
The  difference  in  temperatures  also  means  that  any  electrical  or  thermal  conductivities  or  viscosity 
coefficients  computed  on  the  basis  of  some  combined  temperature  will  be  inaccurate. 

A  two  fluid  model  is  also  neccessary  for  evaluating  the  importance  and  effect  of  plasma  instabilities, 
such  as  those  discussed  by  Choueri  [3],  [4]  and  Hastings  [6] .  These  instabilities  are  often  dependent  on 
the  individual  species  temperatures. 

4.3  Effect  of  Viscosity 

Three  more  cases  were  then  run,  adding  one  effect  in  each  case.  The  first  addition  that  was  made  was 
ambipolar  diffusion.  The  results  for  this  case  are  labeled  Case  2.  The  addition  of  ambipolar  diffusion 
has  caused  a  sharp  decrease  in  the  ionization  fraction,  and  a  smaller  decrease  in  the  gas  temperature. 
The  electron  temperature  has  increased  to  compensate  for  the  diffusion  loss  in  the  ionization  fraction. 
However,  ambipolar  diffusion  has  little  effect  on  the  thruster  performance  as  given  in  Table  5  .  This 
seems  to  be  because  the  electrons  and  ions  which  are  now  being  lost  to  the  side  walls  were  lost  at  the 
channel  exit  in  Case  1.  Although  the  electrodes  will  be  heated  by  the  diffusing  particles,  electrode 
temperature  is  not  included  in  the  model.  However,  as  shown  in  Table  3  for  Case  4,  ambipolar  loss 
absorbs  a  substantial  fraction  of  the  dissipation,  so  that  there  will  be  a  great  deal  of  electrode  heating. 

Recent  work  by  Kilfoyle  [lOj  and  Kuriki  [11]  both  show  measured  gas  temperatures  which  equal  and 
even  exceed  the  electron  temperatures.  Kilfoyle  finds  Tg  ranging  from  leV  up  to  7eV,  or  up  to  80,000K. 
Kuriki  finds  heavy  particle  temperatures  of  up  to  "7,000k.  As  shown  in  Figure  3,  the  inviscid  model 
predicts  gas  temperatures  of  only  5000K.  Some  additional  source  must  be  adding  energy  to  the  heavy 
species.  Heimerdinger  et  al.  [8]  and  Kilfoyle  [10|  propose  viscous  dissipation  as  this  source.  Therefore, 
the  second  effect  examined  was  the  addition  of  viscosity  to  the  model.  The  results  using  this  model  are 
labeled  Case  3.  As  can  be  seen  in  these  figures,  the  heavy  species  temperature  increases  to  the  levels 
found  experimentally.  This  would  seem  to  justify  the  hypothesis  that  viscous  effects  could  cause  high 
heavy  species  temperatures.  The  higher  temperature  leads  to  higher  pressure  throughout  the  channel. 
Also,  the  flow  now  seems  to  be  frictionally  choked,  so  that  the  velocity  levels  off  toward  the  beginning 
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Figure  4:  Two  Fluid  Results:  Velocity 


8 


Z(cm) 

ionization 

recombination 

ambipolar 

diffusion 

0.5 

1.11 

0.26 

0.91 

2.5 

0.68 

0.09 

0.80 

4.5 

1.08 

0.08 

0.78 

6.5 

1.71 

0.10 

0.85 

8.5 

2.25 

0.12 

1.01 

10.5 

2.83 

0.16 

1.22 

12.5 

3.56 

0.20 

1.50 

14.5 

4.58 

0.25 

1.86 

16.5 

6.30 

0.34 

2.37 

18.5 

9.93 

0.56 

3.17 

Table  2:  Magnitude  of  Terms  in  the  Ionization  Fraction  Equation 


of  the  channel,  and  does  not  increase  to  the  levels  found  in  Cases  1  and  2.  As  shown  in  Table  5  the 
thrust  and  the  efficiency  have  both  dropped  considerably,  . 

Heimerdinger  (7j  also  included  viscosity,  but  assumed  thermal  equilibrium.  He  also  computed  the 
ionization  fraction  by  balancing  recombination  with  local  ionization.  His  results  show  an  ionization 
fraction  which  varies  almost  linearly  with  z.  This  is  in  part  due  to  the  large  variation  in  the  overall 
temperature  which  is  used  to  compute  the  ionization.  By  separating  the  electron  and  heavy  species 
temperatures,  this  research  finds  a  lower  ionization  fraction  in  the  bulk  of  the  channel  with  the  ioniza¬ 
tion  fraction  increasing  sharply  at  the  exit.  The  variation  of  the  viscosity  coefficient  and  the  velocity 
distribution  are  also  quite  different  than  in  Heimerdinger. 

One  other  process  which  was  added  to  the  model  was  electron  axial  heat  conduction.  The  results  for 
this  case  are  labeled  Case  4.  The  results  for  this  case  are  quite  similar  to  those  of  Case  3.  The  thrust 
and  efficiency  for  this  case  are  again  shown  in  Table  5. 


4.4  Relative  Importance  of  Effects 

In  steady  state,  the  ionization  fraction  equation  becomes 


da 

dz 


—  miRn,.Snn 
m 


A  c  s  A 

- r  rn,hn, - 

m  rn 


12  D,xnc 

H 2 


(20) 


For  Case  4,  the  complete  model,  the  relative  magnitude  of  each  of  these  effects  at  various  locations 
in  the  channel  is  given  in  Table  2.  It  is  seen  that  recombination  is  a  relatively  weak  effect.  At  the 
beginning  of  the  channel,  ambipolar  diffusion  and  ionization  are  almost  equal.  However,  as  the  electron 
temperature  increases,  ionization  becomes  almost  twice  as  large  as  ambipolar  diffusion.  The  electron 
energy  equation  can  also  be  broken  down  in  this  way.  Again,  in  steady  state, 


m,  J 2 

kpU  a 


m, 

kpU 


E, 


Ei  da  3^,  da  3  dTe  aTe  dp 
k  dz  2  *  dz  2  dz  p  dz 


kpUdz^K'dz'  kPuE'  DaH 2 


In  words,  Ionization  Energy  +  Temperature  Energy  +  Electron  Heating  -  Pressure  Work  =  Dissipation 
-  Collisional  TVansfer  -  Heat  Conduction  -  Ambipolar  Loss.  The  relative  magnitude  of  each  of  these 
terms  is  given  in  Table  3.  For  the  most  part,  dissipation  is  balanced  by  ionization  energy  and  ambipolar 
loss,  although  collisional  energy  transfer  is  also  a  significant  term,  particularly  near  the  inlet.  Notice 
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Dissi¬ 

Collisional 

Heat 

Ambipolar 

pation 

Transfer 

Cond. 

Loss 

3.06 

1.37 

-0.37 

1.67 

1.60 

0.52 

-0.16 

1.46 

2.12 

0.22 

-0,26 

1.41 

3.02 

-0.09 

-0.16 

1.56 

3.64 

-0.48 

-0.15 

1.84 

4.27 

-0.96 

-0.17 

2.24 

5.11 

-1.47 

-0.21 

2.75 

6.53 

-2.0 

-0,28 

3.42 

9.21 

-2  56 

-0.39 

4.34 

15.04 

-3.40 

-0.51 

5.81 

Table  3:  Magnitude  of  Terms  in  the  Electron  Energy  Equation,  xlO  5 

also  that  near  the  exit  the  ions  are  actually  significantly  heating  the  electrons.  Heat  conduction  and 
electron  heating  play  a  small  role,  while  pressure  work  is  negligible. 

5  Area  Variation 

Using  the  full  model  of  Case  4,  the  flow  in  the  constant  area  channel  was  compared  to  the  flow  in 
a  converging-diverging,  or  fully  flared,  channel  (FFC).  The  interelectrode  separation  for  the  channel  is 
shown  in  Figure  5.  The  flow  is  compared  to  that  in  a  constant  area  channel  (CAC)  in  Figures  6-8. 

The  channel  examined  were  chosen  to  be  similar  to  those  examined  experimentaly  by  Heimerdinger, 
Kilfoyle,  and  Martinez  [7],  and  have  the  same  length,  the  same  throat  area,  and  the  same  inlet  and 
exit  area  ratios.  Heimerdinger  proposes  area  variation  as  a  means  of  reducing  inlet  and  exit  current 
concentrations.  This  is  seen  to  be  the  case,  particularly  at  the  exit. 

As  given  in  Table  5,  there  is  a  significant  increase  in  performance  for  the  variable  area  channel. 
This  is  in  part  due  to  a  decrease  in  the  effect  of  viscosity  on  the  velocity,  as  the  viscous  source  term  in 
the  momentum  equation  varies  as  Since  H  has  increased,  the  magnitude  of  this  term  has  dropped, 
despite  the  increase  in  U.  One  other  difference  with  the  constant  area  channel  is  the  higher  electron 
temperature  in  the  fully  flared  channel.  This  can  be  explained  by  examining  Table  4  which  gives  the 
magnitude  of  the  terms  in  the  elect' on  energy  equation.  Comparing  this  table  to  Table  3  points  out 
a  number  of  differences  between  the  two  channels.  First,  collisional  energy  transfer  at  the  ends  of  the 
FFC  is  around  50%  of  its  value  at  the  ends  of  the  CAC.  This  is  because  collisional  energy  transfer  scales 
with  density,  which  is  smaller  at  the  ends  of  the  FFC,  because  of  the  increased  area.  Ambipolar  loss  is 
also  smaller  at  the  ends  of  the  FFC,  because  of  the  increased  interelectrode  distance.  In  the  center  of 
the  channel,  where  both  of  these  effects  are  of  the  same  magnitude  as  in  the  CAC,  the  increased  throat 
dissipation  of  the  FFC  keeps  the  electron  temperature  higher  than  in  the  CAC. 

6  Performance 


The  thrust  and  efficiency  for  the  various  cases  are  shown  in  Table  5.  There  are  a  number  of  loss 
mechanisms  that  prevent  the  thruster  from  being  100%  efficient.  One  reason  for  the  decreased  efficiency 
is  the  frozen  flow  losses.  Since  the  residence  time  of  the  fluid  in  the  channel  is  short,  the  composition  of 
the  fluid  is  not  the  equilibrium  composition.  Therefore,  some  of  the  energy  used  to  ionize  the  neutral 
fluid  is  not  recovered.  Ambipolar  diffusion,  although  it  changes  where  losses  occur,  does  not  seem  to 


Z(cm) 

Ionization 

Energy 

Temper. 

Energy 

Electron 

Heating 

Pressure 

Work 

Dissi¬ 

pation 

Coilisional 

Transfer 

Heat 

Cond. 

Ambipolar 

Loss 

0.5 

2.48 

0.22 

0.0 

-0.17 

4.32 

0.50 

-0.09 

0.94 

2.5 

0.19 

0.02 

0.01 

-0.02 

1.34 

0.23 

-0.09 

0.96 

4.5 

0.78 

0.07 

0.03 

0.01 

0.16 

-0.46 

0.93 

6.5 

2.70 

0.25 

0.02 

0.02 

0.10 

-0.33 

1.17 

8.5 

4.82 

0.46 

0.01 

-0.03 

6.97 

-0.22 

-0.14 

1.83 

10.5 

5.56 

0.53 

0.03 

-0.10 

8.31 

-0.83 

-0.22 

2.74 

12.5 

5.24 

0.51 

0.04 

-0.18 

8.43 

-1.29 

-0.28 

3.48 

14.5 

4.75 

0.48 

0.07 

-0.25 

8.24 

-1.43 

-0.39 

3.81 

16.5 

5.04 

0.52 

0.12 

-0.32 

8.81 

-1.31 

-0.77 

3.84 

18.5 

7.83 

0.85 

0.28 

-0.36 

12.15 

-1.10 

-1.97 

3.83 

Table  4:  Magnitude  of  Terms  in  the  Electron  Energy  Equation,  xlO  5 


Case 

Bo 

ambipolar 

diffusion 

viscosity 

heat 

conduction 

channel 

Et 

V/m 

m 

■ 

1-fluid 

0.1 

no 

no 

no 

CAC,20cm 

433 

4630 

0.62 

2-fluid 

0.1 

no 

no 

no 

CAC,20cm 

399 

4256 

0.57 

2-fluid 

0.1 

yes 

no 

no 

CAC,20cm 

401 

4180 

0.55 

2-fluid 

0.1 

yes 

yes 

no 

CAC,20cm 

258 

2100 

0.21 

2-fluid 

0.1 

yes 

yes 

yes 

CAC,20cm 

287 

2300 

0.23 

2-fluid 

0.1 

yes 

yes 

yes 

CAC,9cm 

469 

3535 

0.34 

2-fluid 

0.1 

yes 

yes 

yes 

FFC,9cm 

600 

4717 

0.47 

Table  5:  Thrust  and  Efficiency  for  All  Cases 


change  the  amount  of  loss.  However,  adding  viscosity  to  the  model  does  cause  a  large  drop  in  efficiency. 
At  the  low  inlet  magnetic  field,  heat  conduction  does  not  affect  the  efficiency  significantly. 

7  Conclusions 

This  research  had  a  number  of  goals.  The  first  was  to  develop  a  numerical  method  with  which 
to  analyze  a  two  fluid  model.  That  goal  has  been  achieved,  and  it  has  been  shown  that  there  are 
some  significant  differences  between  the  one  fluid  and  the  two  fluid  results.  The  method  developed  has 
proven  to  be  very  versatile,  and  should  become  a  useful  tool  in  thruster  analysis.  Also,  the  method  has 
been  tested  against  previous  one  fluid  work  and  has  been  shown  to  be  accurate.  A  second  goal  was  to 
determine  the  effect  of  viscosity  on  the  thruster  fluid  flow.  It  has  been  seen  that  viscosity  raises  the 
gas  temperature  to  levels  fourd  experimentally.  The  relative  importance  of  various  effects  on  the  flow 
ha3  also  been  calculated.  Finally,  area  variation  has  been  shown  to  increase  efficiency  by  a  significant 
amount,  and  reduce  exit  current  concentrations. 
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Abstract 

A  two-dimensional  numerical  model  has  been  developed 
in  order  to  analyze  electro-magnetic  plasma  accelerators 
also  called  Self-Field  Magneto-Plasma-Dynamic  Thrusters. 
This  model  uses  a  Magneto-Hydro-Dynamic  description  of 
the  gas  considered  as  a  fully  ionized,  isothermal  plasma, 
and  takes  into  account  the  Hall  effect  (non  linear  conduc¬ 
tivity)  and  the  interaction  between  the  magnetic  field  and 
the  fiuid  dynamics  of  the  plasma.  The  system  of  equations 
is  discretized  into  finite  volumes,  and  is  solved  by  a  Newton- 
Raphson  scheme.  Results  from  the  MHD  model  were  calcu¬ 
lated  for  a  mass  flow  rate  of  6  g/s  of  Argon  and  for  currents 
up  to  ten  kilo- Amperes. 

Introduction 

During  the  past  years  the  development  of  high  current 
plasma  accelerators  (also  called  Magneto-Plasma- Dynamic 
Thrusters  )  has  been  aimed,  among  other  things,  at  im¬ 
proving  the  efficiency  of  energy  conversion  from  electrical 
to  kinetic  energy.  Several  researchers  have  observed  that 
the  efficiency  increases  with  the  current,  but  at  some  value 
a  transition  occurs  to  an  unstable  and  destructive  regime  of 
the  discharge  (called  onset  of  instability  in  the  literature), 
when  the  ratio  exceeds  a  critical  value  which  depends 
on  the  nature  of  the  propellant  and  the  thruster  geome¬ 
try.  Some  researchers  [3)  (llj  have  proposed  a  theory  for 
the  transition  to  an  unstable  behavior  based  on  the  deple¬ 
tion  of  the  anode  from  current-carrying  electrons,  the  for¬ 
mation  of  a  starvation  layer  along  the  anode  characterized 
by  current  densities  in  excess  of  the  random  thermal  flux  of 
electrons,  and  the  reversing  of  the  potential  drop  between 
the  the  plasma  and  the  electrode.  In  order  to  study  the 
formation  of  this  starvation  layer  one  needs  to  take  the  geo¬ 
metrical  properties  of  the  discharge  into  account,  or  in  other 
words,  to  look  at  a  two-dimensional  model  of  the  discharge. 
Eventually  this  model  can  be  used  to  find  optimal  shapes 
that  maximise  certain  parameters  like  efficiency  or  electrode 
life. 
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Various  attempts  have  been  made  to  develop  numerical 
two-dimensional  models  of  MPD  discharges  using  finite  el¬ 
ement  or  finite  difference  descriptions.  I.  Kimura  et  a l.  [17] 
have  presented  a  two  dimensional  model  for  the  electromag¬ 
netic  field  (self  field),  including  the  Hall  effect.  The  model 
was  primarily  designed  to  calculate  the  current  distribution 
near  the  cathode.  They  later  extended  their  model  [27]  to 
systems  that  work  with  an  external  focussing  magnetic  field. 
These  models,  however,  were  not  designed  for  an  accurate 
modelling  of  realistic  thruster  geometries,  and  did  not  in¬ 
clude  the  effect  of  the  convection  of  the  electromagnetic  field 
by  the  plasma.  This  last  limitation  was  partially  removed 
in  a  later  paper  [28]  where  a  simplified  calculation  of  the 
Sow  was  made  and  coupled  with  the  magnetic  field.  Ao  and 
Fujiwara  Jl]  have  attempted  to  develop  a  one- fluid  Magneto- 
Hydro-  Dynamic  model  of  the  plasma  where  the  equations 
describing  the  motion  of  the  plasma,  considered  as  a  neu¬ 
tral  conductive  fiuid,  are  coupled  with  the  magnetic  field 
through  the  Lorenti  force  and  the  ohmic  heating.  However 
this  model  was  limited  to  a  very  simple  geometry. 

It  was  decided  to  develop  a  two-dimensional  numerical 
model  based  on  a  MHD  description  of  the  plasma  in  the 
collision  dominated  regime,  in  order  to  calculate  the  steady 
state  distribution  of  the  electro-magnetic  field  and  the  mo¬ 
tion  of  the  plasma  in  the  accelerator.  We  focused  our  atten¬ 
tion  on  the  phenomena  observed  at  the  electrodes,  especially 
at  the  anode,  in  an  effort  to  detect  the  onset  of  the  turbulent 
regime. 

The  system  of  partial  differential  equations  is  discretized 
using  a  finite  volume  technique  used  extensively  in  Compu¬ 
tational  Fluid  Dynamic  algorithms.  The  system  of  equa¬ 
tions  was  solved  using  a  Newton-Raphson  method,  derived 
from  an  initial  calculation  of  the  magnetic  field  alone  for 
which  this  method  works  well  due  to  the  dominant  elliptic 
behavior  of  the  equations.  The  method  of  discretization  of 
the  equations  (finite  volume  method)  as  well  as  the  idea 
of  the  addition  of  an  artificial  viscosity,  was  derived  from 
Jameson  and  Yoon  [15].  A  confirmation  of  the  possibility 
of  the  Newton-Raphson  method  was  found  in  the  lecture  of 
Jeepenen  [16],  as  well  as  in  the  doctoral  thetee  of  M.  Giles 
[8]  and  M.  Drela  [7| 
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Physical  Model 


The  equation*  deecribing  the  model  are  derived  from  the 
general  Magneto-Hydro-Dynamic  approximation  which  in¬ 
clude*  the  following  equation*:  conservation  of  mass,  mo¬ 
mentum  and  energy  in  the  inviscid  fluid;  the  Ohm’s  law 
(derived  from  the  electron  momentum  conservation);  and 
the  quasi-static  Maxwell’s  equations. 

3,p  +  V(pu)  =  0 
3,pn  -t-  Div(^uu  -  S  -  M)  -  0 

Where  5  is  the  pressure  stress  tensor,  and  ,Vf  is  the  Maxwell 
stress  tensor. 

3,{EeK<  ^  E,  ~  EK)  +  V  {[E,  +  EK  +  P)  =  0 

fJo 

Where  Eem,  E/,  Ek  are  respectively  the  electromagnetic 
energy  (  B3/ 2^io),  the  internal  energy,  and  the  kinetic  en- 

p  D 

ergy,  per  unit  volume;  F  is  the  pressure;  ^  is  the  Poynt- 
ing  vector;  q  is  the  heat  conduction.  The  energy-equation  is 
usually  completed  by  a  state  equation  relating  the  pressure 
to  the  internal  energy,  for  instance  : 

P  =  h- 1  )E, 

The  generalised  Ohm’s  law  can  be  written  as: 

1  1  VP 

E+uxB=-J  +  — -J  x  B - — 

cr  en,  en. 

The  Maxwell’s  equations  reduce  to: 


V  x  E  =  -3,  B  VB  =  0 

V  x  B  =  pqJ 


We  have  assumed  that  the  field  is  quasi-static  and  that  the 
plasma  is  approximately  neutral. 

Let  us  now  look  at  the  possible  simplifications: 

Geometry:  The  flow  and  the  electromagnetic  field  are 
axisymetric.  There  is  no  externally  imposed  magnetic  field. 
The  Magnetic  field  is  therefore  azimuthal. 


The  Ohm’s  Law  and  the  electromagnetic  equations  are 
combined  in  order  to  eliminate  E  and  J.  One  obtains  an 
equation  in  B  that  describes  the  diffusion  and  convection  of 
the  magnetic  field  in  a  conducting  medium: 


_  ,  _  V  x  B  (V  x  B)  x  B  VP., 

V  x  (-u  x  B  + - +  i - - - )  =  -3, 


Po  en. 


Rewriting  the  equation  in  non-dimensional  form, 


-d.’B*  +  V*  x  (n*  x  B*)  = 

1  .  fV’xB-  (VxB-)xB’  /?V(n;77) 

Fin  (  <7-  n;  2  n; 


introduce*  three  non-dimensional  number*:  M  =  -J—  for 
the  magnetic  diffusion  due  to  conductivity;  Ni  =  #*-  for 
the  non  Linear  Hall  effect;  Ns  =  for  the  electron  dif¬ 
fusion.  ftn  =  fio&oVoLo  is  the  magnetic  Reynolds  number, 
based  on  the  distance  between  the  electrodes;  H,  is  the  Hall 
parameter; 


H, 


Bo<To 

en,o 


wc.r 


where  <jr,  is  the  electron  cyclotron  angular  frequency  and 
r,  =  l/i/,  is  the  electron  collision  time  .  0  is  the  plasma- ,3; 
ratio  of  the  pressure  to  the  magnetic  pressure, 


0  = 


2po  P. 

Bl 


The  values  of  these  non-dimensional  numbers  give  an  indica¬ 
tion  of  the  relative  importance  of  each  phenomenon:  Inside 
the  channel  Ni  s»  0.5,  Ni  «  0.6,  Nj  m  0.03.  Therefore  we 
take  the  Hall  effect  and  the  convection  into  account,  but  we 
neglect  the  effect  of  the  electron  pressure  gradient  on  the 
conductivity.  (Notice  that  a  magnetic  Reynolds  Number 
based  on  the  channel  length  would  be  of  order  5  to  10). 


Constant  temperature  approximation:  In  a  plasma  like 
the  one  considered  here,  the  energy  equation  is  the  most 
difficult  to  justify.  At  high  temperature  (around  12500 
K)  the  gas  is  strongly  ionized  and  radiates  energy.  These 
processes  are  not  necessarily  in  equilibrium  and  therefore 
cannot  be  accurately  described  by  the  energy  conservation 
equation  used  in  the  fluid  theory  at  low  temperatures.  As  in 
most  MHD  models  the  energy  conservation  equation  will  be 
dropped  and  replaced  with  the  assumption  that  the  plasma 
is  roughly  isothermal.  This  assumption  can  be  justified  by 
the  effect  of  ionization  which  tends  to  oppose  any  change 
in  temperature  by  a  displacement  of  the  ionization  equilib¬ 
rium,  by  the  thermal  conductivity  of  the  plasma,  and  by 
the  value  of  the  pressure  relative  to  the  magnetic  pressure 
0  <C  l.  The  non-dimensional  number  describing  the  balance 
between  kinetic  energy  and  heat  conduction  (equivalent  of 
the  Reynolds  number  for  the  energy  equation)  is: 


Ci 


LqPqUq 

xTo 


Its  value  goes  from  as  200  :n  the  cathode  jet,  to  as  1  along 
the  external  surface  of  the  anode. 


Conductivity:  In  the  Spitzer-Harm  aproximation  the  con¬ 
ductivity  of  a  fully  ionized  plasma  is  proportional  to  r*. 
Consequently  we  have  a  constant  value  for  the  electrical 
conductivity  a  —  3800Si/m. 

Even  though  the  enrgy  balanc«  is  not  written,  the  fact 
that  one  part  of  the  electrical  work  E  J  is  dissipated  ohmi- 
cally  as  is  implicit  in  the  rest  of  the  formulation,  and 
therefore  a  meaningful  thruster  efficiency  can  be  computed. 
What  is  ignored  is  the  ultimate  fate  of  the  ohmic  heat  and 
it*  apportionment  to  ionisation,  heat  conduction,  radiation, 
etc. 
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Since  ionization  processes  are  not  covered  in  the  model, 
the  plasma  ie  comidered  fully  ionised  when  it  enters  the 
inlet  of  the  accelerator. 

Although  the  read  Reynolds  number  is  of  the  order  of  100 
to  200,  the  plasma  is  considered  inviscid,  which  means  that 
second  order  derivatives  appear  only  in  the  electromagnetic 
equation. 

The  final  set  of  equations  after  simplifications  is  then; 
3,(pu.)  +  id,(rpu,)  =  0 

Q 2  1 

d,(pv,u,  -t -  P  +  - — )  +  -dr(rpu.Ur)  =  0 

2/io  r  ' 

3,{p U.Ur)  -r  -dr(r(pUrUr  +  P  +  - - ))  +  - ^2.  =  0 

r  2uo  r 

P  =  pRTo 

Edl  =  0 

E.  =  —  -3rrB - ~d,B  -  urB 

MoCT  r 

Er  =  -  —  a.B - —  -d,rB  +  u  ,B 

Ho<7  Vkitn.  r 

The  boundary  conditions  were  taken  as  follows: 


The  existence  of  a  steep  density  gradient  immediately  ad¬ 
jacent  to  the  back  plate  makes  it  very  difficult  to  resolve  this 
layer  numerically.  Fortunately  one  of  us  has  shown  (22)  that 
the  details  of  the  inlet  gradients  have  only  a  small  influence 
downstream.  In  our  calculation  we  assign  an  arbitrary  value 
to  p. 

On  the  outer  boundary  of  the  computational  domain  (i.e. 
far  into  the  plume)  B  =  0.  This  condition  means  that  there 
is  no  current  going  outside  the  boundary  of  the  calculation, 
which  is  therefore  an  insulating  surface.  There  are  no  gra¬ 
dients  of  p,  pu.,  p ur  along  the  direction  parallel  to  the  flow. 

On  the  axis  of  symetry  B  =  0. 

Numerical  Method 

The  system  of  equations  is  solved  by  the  N’ewton- 
Raphston  method  using  the  technique  of  Finite  Volumes 
for  the  discretization  of  the  equations  of  motion,  and  a  tech¬ 
nique  similar  to  the  method  of  finite  volumes  for  the  mag¬ 
netic  field.  Artificial  viscosity  is  added  to  the  equations  of 
motion  in  order  to  avoid  non-linear  instabilities  at  places 
of  strong  density  gradients.  The  conservation  equations  are 
discretized  according  to  the  method  of  finite  volumes:  The 
computational  space  is  divided  into  20  x  41  quadrilaterals. 
The  grid  is  shown  in  figure  1.  The  variables  (p,pu,,pu,,B) 
are  approximated  by  their  space-average  on  each  quadrilat¬ 
eral.  The  equations  are  written  in  their  integral  form: 


The  tangential  electric  field  is  zero  on  the  electrodes.  The 
electrodes  include  the  conducting  surfaces  of  the  anode  and 
the  cathode.  The  voltage  drop  across  the  plasma-sheath, 
which  extends  over  a  few  Debye  lengths  inside  the  plasma, 
will  be  neglected.  Consequentely  the  voltage  drop  across  the 
channel  will  be  underestimated  when  the  anode  is  starved, 
i.e.  at  high  current  densities. 

The  entrance  surface  (the  back  plate)  is  an  insulating 
surface  which  carries  the  injector  holes.  It  is  assumed  that 
the  gas  is  injected  uniformly  across  the  back  plate,  which 
would  be  the  case  with  a  porous  back  plate  or  with  a  plate 
containing  a  large  number  of  injector  holes.  A  transition 
between  the  gas  state,  characterized  by  a  low  temperature 
and  a  low  conductivity,  and  the  plasma  state,  where  the 
gas  is  ionized  and  is  a  good  conductor  of  electricity,  arises 
once  the  gas  is  in  the  annular  chamber.  This  ionization 
takes  place  in  a  region  close  to  the  back  plate  surface.  The 
calculation  starts  after  the  ionization  layer  assuming  that 
the  gas  has  reached  a  uniform  temperature  of  12500  K,  at 
which  the  electrical  conduction  is  sufficently  high  to  consider 
the  gas  as  a  plasma.  Across  the  entrance  plane  the  magnetic 
field  is:  B  =  noI/2wr.  The  dynamic  variables  p,  pu,,puf 
are  specified  along  the  entrance  surface  :  pu,  =  0;  pu,  is 
specified  by  the  mass  flow  rate: 

m 

PM.  =  - - 


where  the  the  tensor  puti  —  S  —  A/  has  the  following  physical 
components  in  cylindrical  coordinates  (z,r,  <f>)  : 


Since  the  problem  is  axisymetrical,  the  volume  elements 
are  obtained  by  rotating  the  quadrilaterals  that  define  the 
grid  cells  by  an  infinitesimal  angle  Ad  around  the  axis  of 
symetry.  The  surface  integral  around  each  element  is  bro¬ 
ken  into  six  integrals,  one  for  each  facet.  Each  of  them 
is  discretized  by  assuming  that  the  fluxes  are  constant  on 
each  facet  of  the  boundary.  For  the  electro-magnetic  field 
one  uses  a  different  integral  (line  integral  over  a  closed  con¬ 


tour): 


£ 


Edl  =  0 


The  contours  are  defined  by  the  same  quadrilaterals  as 
above.  For  each  iteration  the  Newton-Raphson  algorithm 
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invert*  a  matrix  of  order  3200  which  represent*  the  *y*tem 
of  partial  differential  equation*  Unearned  at  the  current  e»- 
timate  of  the  solution  vector.  Thi*  take*  about  10  minute* 
on  a  microVax.  One  need*  about  10  iteration*  to  attain 
convergence. 

Experimental  Background 

Experimental  studies  have  identified  three  modes  of  op¬ 
eration  for  the  plasma  accelerator:  electrothermal  mode, 
for  which  the  dominant  acceleration  mechanism  is  the  ex¬ 
pansion  of  the  heated  gas  across  the  thermodynamic  pres¬ 
sure  gradient;  electromagnetic  mode,  for  which  the  dom¬ 
inant  acceleration  mechanism  is  the  expansion  of  the  gas 
across  the  magnetic  pressure  gradient;  finally  the  unstable 
regime  which  set  on  when  the  critical  ratio  is  exceeded. 
For  argon  as  propellant,  Hugel  (ll|  give*  an  empirical  value 
for  the  limiting  current, 

!—  ss  2.5  10l°  (A3  kg-1  s)  , 

m 

which  gives  a  current  of  the  order  of  12  kA  for  rh  =  6  g/s. 
More  recent  work  at  Princeton  and  elsewhere  ha*  raised 
this  by  factors  of  up  to  3,  depending  on  the  geometry  and 
the  injector  configuration.  An  onset  prediction  based  on 
the  total  starvation  of  the  anode  was  proposed  by  Baksht, 
Moishes  and  Rybakov  (3|, 

fm]  _  945  tkT  LHa 
rh*  8  crfi*  d& 

which  leads  to  a  critical  value  of  19  kA  for  our  geometry. 

An  explanation  which  does  not  include  the  Hall  effect 
has  been  proposed  by  Lawless  and  Subramanian  [19],  based 
on  the  mechanism  of  an  excessive  back  EMF  induced  by  a 
strong  convective  effect  a  x  B  in  the  opposite  direction  of 
the  appUed  driving  field  £.  Neglecting  the  Hall  effect  they 
obtained,  from  a  one-dimemsional  analysis,  an  expression 
for  the  limiting  current: 

—r*  as  8.52—  (A3  kg-1  «) 

m  /io 

where  a*  is  the  speed  of  sound  at  the  choking  point  and  jio  is 
the  permeabiUty  of  the  vacuum.  This  give*  a  current  of  the 
order  of  20  kA  for  a*  =  104  m/s  and  m  =  6  g/s.  However 
the  physical  basis  of  their  argument  isopen  to  criticism  [22]. 

Assuming  that  the  magnetic  field  is  confined  between  two 
coaxial  cylindrical  electrodes  and  doe*  not  extend  outside 
the  exit  plane  of  the  thruster,  one  can  Integrate  the  effects 
of  the  Maxwell  stress  tensor,  and  one  can  predict  the  thrust 
due  to  the  electro-magnetic  effects  by  the  formula: 

T.m  =  ~-  In  - 

4*  r, 

(11  Newton  at  lOkA  for  ^  =  3)  The  ratio  give*  the 
exit  velocity  e,. 


The  voltage  is  expected  to  be  linear  with  the  current  for 
low  current*  (Joule  heating),  and  to  increase  at  a  higher 
rate  for  high  currents  as  the  field  accelerate*  the  plasma 
and  competes  against  the  induced  back  EMF,  showing  a 
cubic  dependence  with  the  intensity. 

Results 

We  have  run  fully  coupled  calculations  for  a  mass  flow  rate 
of  6  grams  per  second  and  for  various  electrical  currents  be¬ 
tween  0  and  lOkA.  Uncoupled  calculations  (i.e.  calculating 
the  electromagnetic  field  with  a  frozen  flow  field)  have  been 
obtained  for  currents  between  0  and  100  kA. 

A  maximum  of  the  magnetic  field  is  observed  along  the 
cathode.  Figure  2  shows  the  contours  of  constant  rB.  For  an 
axisymetrical  geometry  these  lines  also  represent  the  lines 
of  current.  The  influence  of  the  magnetic  field  on  the  con¬ 
ductivity  tensor  (Hall  effect),  is  seen  in  the  way  the  current 
lines  concentrate  at  the  tip  of  the  anode  and  at  the  root 
of  the  cathode.  As  the  current  increases  this  effect  becomes 
more  pronounced.  This  effect  can  also  be  observed  on  figure 
3,  which  shows  the  current  density  at  the  electrodes. 

Figure*  4  and  5  show  the  density  contours  for  0  and  lOkA. 
For  lOkA  (figure  5)  the  density  shows  strong  gradients  and 
large  fluctuations.  At  the  inlet  the  density  drops  abruptly 
over  a  few  cells.  The  compression  of  the  plasma,  under  the 
influence  of  the  Lorents  force  (inward  pinching  effect),  can 
be  observed  along  the  cathode  and  along  the  axis  of  symetry. 
Figure  0  shows  an  enlargement  of  the  density  plot  for  lOkA, 
where  one  observes:  the  compression  of  the  plasma  along  the 
cathode;  the  strong  rarefaction  that  follows  the  expansion 
of  the  plasma  around  the  cathode  corner,  which  reaches  a 
minimum  of  1.5  IO30  m~3;  the  recompression  along  the  axis 
of  symetry,  where  the  density  reaches  a  maximum  of  108. 
IO30  m“3  at  a  distance  of  1  cm  from  the  tip  of  the  cathode. 
On  the  anode  side  the  plasma  expands  monotonically  from 
33.  IO30  m~3  to  0.15  IO30  m~3 .  Figure  7  shows  the  density 
across  the  thruster  along  three  different  rows  of  cells.  It 
shows  the  considerable  increase  in  density  in  the  cathode 
jet,  downstream  of  the  cathode  along  the  axis  of  symetry. 

A  plot  of  the  velocity  vector  obtained  for  lOkA  (figure 
8  )  shows  a  strong  acceleration  in  the  surface  layer  that 
surrounds  the  cathode.  The  existence  of  this  layer  can  be 
inferred  from  the  high  values  of  both  the  electric  field  and 
the  current  density  along  the  cathode  (about  an  order  of 
magnitude  higher  than  at  the  anode). 

The  plot  of  the  Hall  number  ■—  (figure  9)  shows  con¬ 
siderable  gradient*  on  the  anode  and  cathode  tips.  These 
gradients  are  related  to  the  electron  density  gradients  cre¬ 
ated  by  the  expansion  of  the  plasma  as  it  follows  the  elec¬ 
trode  surface*.  The  Hall  effect  (Tensorial  conductivity)  is 
likely  to  be  responsible  for  the  computational  instabilities 
encountered  at  higher  current*  along  the  electrode  surfaces. 
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The  thrust  was  computed  for  various  currents  and  is  listed 
in  table  2.  The  thrust  increase  due  to  the  current  is  appar¬ 
ently  very  modest  (9.8  N  at  lOkA),  but  is  comparable  to 
the  thrust  due  to  the  Maxwell  stress  tensor  predicted  by 
the  formula  : 


(11  N  at  10kA  for  =  3  ).  Thus  the  thruster  is  still  oper¬ 
ating  at  this  condition  mainly  in  the  electrothermal  regime. 

Conclusions 

The  calculation  of  the  MHD  model  was  done  for  three  val¬ 
ues  of  the  current:  0,  5kA  and  lOkA.  The  model  does  not 
converge  above  lOkA.  The  o-obable  cause  of  the  computa¬ 
tional  instability  is  the  exct  ve  value  of  the  Hall  parameter 
along  the  surfaces  where  rarefaction  is  present.  Above  a  cer¬ 
tain  value  of  the  current  this  rarefaction  will  be  increased 
by  the  action  of  the  Lorent*  force,  leading  to  a  failure  of  th' 
model  along  the  boundaries. 

The  same  cause  may  be  responsible  for  the  onset  of  insta¬ 
bilities  observed  in  experimental  devices.  (The  onset  cur¬ 
rent  expected  from  Hugel’s  formula  is  12k A,  the  value  from 
the  Baksht  et  al.  is  19kA,  and  the  value  from  Lawless  and 
Subramanian  is  20kA). 

The  artificial  damping  introduced  in  the  equation  of  mo¬ 
tion  in  order  to  obtain  a  reasonably  well-posed  problem 
modifies  the  solution  by  an  amount  that  reaches  5  percent 
of  the  momentum  flux  (for  lOkA).  The  largest  contribution 
to  this  error  appears  at  the  sharp  turn  around  the  anode 
corners. 
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Symbols 


a"  Speed  of  sound  at  chocking  point  in  the 

theory  of  Lawless  and  Subramanian 
B  Magnetic  field  (T) 

c.  Exit  velocity 

Ci  Non  dimensional  parameter 

in  the  energy  equation  (Clausius  Number) 
i  distance  between  electrodes 

t  Charge  of  the  proton 

E  Electric  field  (  V  m~l) 

Eem  Electromagnetic  energy  per  unit  volume 

Ei  Internal  energy  per  unit  volume 

Ek  Kinetic  energy  per  unit  volume 

H  Average  perimeter  of  the  inlet 

Ha  Hail  parameter 

I’  Current  at  the  threshold 

of  the  turbulent  regime. 

J  Current  density  (  A  m~3  ) 

k  Boltzman  constant 

L  Average  length  of  the  channel 

m,  Molecular  mass  of  the  gas  (  kg) 

Af  Maxwell  stress  tensor  (  N  m-*) 

rh  Mass  flow  rate  of  the  propellant  (kg  s-1) 

m„m,  Z-  and  r-momentum  (=  pu,,  pu,) 

n.  Density  of  electrons  (  m-s  ) 

P  Pressure  (  N  m_a  ) 

q  Heat  flux  (  W  m~*  ) 

R  Gas  constant  (  J  kg-1  K~l  ) 

Magnetic  Reynolds  number 
5  Stress  tensor  (  N  m~3) 

T  Temperature  (  K  ) 

u  Velocity  (  m  s-1  ) 

uu  Diad  (u,u,  in  cartesian  coordinates) 


0  Plasma- Beta  (pressure/magnetic  pressure) 

7  Specific  heat  ratio 

k  Heat  conductivity 

Ho  Permittivity  of  vacuum 

p  Mass  density  (  kg  m-3  ) 

tr  Conductivity  (  Si  m-1  ) 

t.  Electron  collision  time 

<t>  Axisymertic  coordinate 

Electron  gyrofrequency  (angular) 
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Cathode  radius 

rc 

- - - 

2.  cm 

Anode  radius 

6.  cm 

Cathode  length 

6.  cm 

Anode  length 

12.  cm 

Inlet  density 

Ptnlet 

2.21  1CT4  kg  m~3 

Inlet  pressure 

Ptnlet 

1144  Pa 

Inlet  velocity 

Z inlet 

2700  m/s 

In. let  X-sec.  area 

WinJet 

cs 

e 

r* 

l 

o 

Temperature 

Tu 

1.25  104  K 

Ctnductivity 

<?Q 

3.8  103  Si/m 

Mass  flow  rate 

•n 

6  g/« 

Current 

(0.,  5.,  10.)  kA 

Gas  constant 

R 

0.416  kJ/kg/K 

Table  1:  Parameters  for  the  Calculations 


RESULTS  OF  THE  MODELS 


Current 

Voltage 

Thrust 

Isp 

' 

MHD  model 

(MHD) 

(MHD) 

!  (kA) 

(V) 

(N) 

(km/s) 

0. 

30.2 

5.0 

j  5. 

6.7 

32.2 

. 

5.4 

10. 

' 

16.2 

40. 

6.7 

Table  2:  Results 
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Figure  2:  Current  lines  for  I  =  10*A  (Calculation) 
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Figure  6:  Density  for  /  =  IQkA  (Calculation) 


